function [S, dSdpar] = T1IRAbsMagnitude( par, echoes,TR)

% computes the magnitude M with given parameters and inversion times
% echoes = [n x 1] inversion times
% par = the 3 element parametervector [A, B, R1]'
%
% The datamodel is given by
% M = A - B*exp(-echoes*R1).
%
% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007

if size(par,1) ~= 3 %HENK
    error('incorrect parameter vector');
end;

A = par(1);
B = par(2);
R1 = par(3);
t = echoes;

% %mckinzey with TR
E1 = exp(min(650,-t*R1));
E2 = exp(min(650, -TR*R1));
S1 = A*(1-B*(E1)+E2);
S = abs(S1);

if nargout>1   
    sgn = sign(S1);
    
%     %mckinzey with TR
    dSdR1 = ((t*A*B) .*E1 - A.*(TR).*E2) .* sgn;
    dSdA  = (1 - B.*E1 + E2) .* sgn;
    dSdB  = -A.*E1.*sgn;
    dSdpar = [dSdA dSdB dSdR1];
end
    
%old implement without McK term    
% numtr = size(par,2);
% numgr = size(echoes,1);
% 
% A = zeros(numgr,numtr);
% A_ex = A;
% A_part = A;
% Aabs = A;
% 
% for k=1:numtr
%  A_ex(:,k) = exp(-echoes(:,1)*par(3,k)); %HENK matrix with all results 
%  A_part(:,k) = par(2,k) * A_ex(:,k);
%  Aabs(:,k) = par(1,k) - par(2,k) * A_ex(:,k);   
%  A(:,k) = abs(par(1,k) - par(2,k) * A_ex(:,k));   
% end;
% 
% if nargout>1
%     dAdpar=([A(:,1) A(:,1) A(:,1)]./[Aabs(:,1) Aabs(:,1) Aabs(:,1)]).*[ones(size(A_ex,1),size(A_ex,2)) -A_ex(:) repmat(echoes,numtr,1).*A_part(:)]; %HENK 
% end;
